library("reshape2")
library("ggplot2")
theme_set(theme_bw(base_size = 12) +
theme(rect = element_rect(fill = "transparent")))
library("leaflet")
library("rgdal")
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.0, released 2018/12/14
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ shared files: (autodetected)
## Linking to sp version:1.4-5
library("htmlwidgets")
library("lavaan")
## This is lavaan 0.6-7
## lavaan is BETA software! Please report any bugs.
library("lavaanPlot")
rm(list = ls())
setwd("~/epi")
These helper functions retrieve, read and process data.
source("UtilityFunctions.R")
Run the helper function importEPI() and release the contents of the resulting list into the global environment.
epi.list <- importEPI()
## 'data.frame': 8234 obs. of 7 variables:
## $ code : int 4 24 8 784 32 51 28 36 40 31 ...
## $ iso : chr "AFG" "AGO" "ALB" "ARE" ...
## $ country : Factor w/ 179 levels "Afghanistan",..: 1 4 2 170 6 7 5 8 9 10 ...
## $ region : Factor w/ 8 levels "Asia-Pacific",..: 7 8 2 5 6 3 6 4 4 3 ...
## $ X2019 : num 1.93e+10 8.88e+10 1.53e+10 4.21e+11 4.45e+11 ...
## $ EPI.new : Factor w/ 46 levels "AGR","AIR","APE",..: 12 12 12 12 12 12 12 12 12 12 ...
## $ EPI.new.value: num 25.5 29.7 49 55.6 52.2 52.3 48.5 74.9 79.6 46.5 ...
list2env(epi.list, .GlobalEnv)
## <environment: R_GlobalEnv>
rm(epi.list)
Download the shape file, if it is absent. Note: overwrites existing files!
if(!file.exists("2020/TM_WORLD_BORDERS-0.3.shp")){
print("File does not exist, downloading it now:")
download.file("thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip",
destfile="world_shape_file.zip")
system("unzip -u world_shape_file.zip;
rm world_shape_file.zip Readme.txt")
}
The, read this shape file into a SpatialPolygonsDataFrame.
world.spdf <- readOGR(dsn = "2020", layer = "TM_WORLD_BORDERS-0.3",
verbose = TRUE)
## OGR data source with driver: ESRI Shapefile
## Source: "/home/r/epi/2020", layer: "TM_WORLD_BORDERS-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings: POP2005
summary(world.spdf)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -180 180.0000
## y -90 83.6236
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
## FIPS ISO2 ISO3 UN
## Length:246 Length:246 Length:246 Min. : 4.0
## Class :character Class :character Class :character 1st Qu.:215.0
## Mode :character Mode :character Mode :character Median :429.0
## Mean :431.8
## 3rd Qu.:650.5
## Max. :894.0
## NAME AREA POP2005 REGION
## Length:246 Min. : 0.0 Length:246 Min. : 0.00
## Class :character 1st Qu.: 44.5 Class :character 1st Qu.: 2.00
## Mode :character Median : 5515.5 Mode :character Median : 19.00
## Mean : 52696.1 Mean : 65.43
## 3rd Qu.: 34708.8 3rd Qu.:142.00
## Max. :1638094.0 Max. :150.00
## SUBREGION LON LAT
## Min. : 0.00 Min. :-178.13 Min. :-80.4460
## 1st Qu.: 14.00 1st Qu.: -50.16 1st Qu.: -0.3025
## Median : 30.00 Median : 17.66 Median : 16.5110
## Mean : 54.84 Mean : 13.29 Mean : 16.4075
## 3rd Qu.: 61.00 3rd Qu.: 50.01 3rd Qu.: 39.1067
## Max. :155.00 Max. : 179.22 Max. : 78.8300
Unfortunately, as soon as it comes to map the world, geopolitics come to the fore. This is reflected in the mismatches between country lists due to the different recognition statuses of these countries.
For one, there are distinctly more countries in the shp file:
length(world.spdf$NAME)
## [1] 246
length(epi$country)
## [1] 179
These mismatches stem mainly from different naming schemes.
sort(setdiff(world.spdf$NAME, epi$country))
## [1] "Åland Islands"
## [2] "American Samoa"
## [3] "Andorra"
## [4] "Anguilla"
## [5] "Antarctica"
## [6] "Aruba"
## [7] "Bermuda"
## [8] "Bouvet Island"
## [9] "British Indian Ocean Territory"
## [10] "British Virgin Islands"
## [11] "Burma"
## [12] "Cape Verde"
## [13] "Cayman Islands"
## [14] "Christmas Island"
## [15] "Cocos (Keeling) Islands"
## [16] "Congo"
## [17] "Cook Islands"
## [18] "Democratic Republic of the Congo"
## [19] "Falkland Islands (Malvinas)"
## [20] "Faroe Islands"
## [21] "French Guiana"
## [22] "French Polynesia"
## [23] "French Southern and Antarctic Lands"
## [24] "Gibraltar"
## [25] "Greenland"
## [26] "Guadeloupe"
## [27] "Guam"
## [28] "Guernsey"
## [29] "Heard Island and McDonald Islands"
## [30] "Holy See (Vatican City)"
## [31] "Hong Kong"
## [32] "Iran (Islamic Republic of)"
## [33] "Isle of Man"
## [34] "Jersey"
## [35] "Korea, Democratic People's Republic of"
## [36] "Korea, Republic of"
## [37] "Lao People's Democratic Republic"
## [38] "Libyan Arab Jamahiriya"
## [39] "Liechtenstein"
## [40] "Macau"
## [41] "Martinique"
## [42] "Mayotte"
## [43] "Micronesia, Federated States of"
## [44] "Monaco"
## [45] "Montserrat"
## [46] "Nauru"
## [47] "Netherlands Antilles"
## [48] "New Caledonia"
## [49] "Niue"
## [50] "Norfolk Island"
## [51] "Northern Mariana Islands"
## [52] "Palau"
## [53] "Palestine"
## [54] "Pitcairn Islands"
## [55] "Puerto Rico"
## [56] "Republic of Moldova"
## [57] "Reunion"
## [58] "Saint Barthelemy"
## [59] "Saint Helena"
## [60] "Saint Kitts and Nevis"
## [61] "Saint Martin"
## [62] "Saint Pierre and Miquelon"
## [63] "San Marino"
## [64] "Somalia"
## [65] "South Georgia South Sandwich Islands"
## [66] "Svalbard"
## [67] "Swaziland"
## [68] "Syrian Arab Republic"
## [69] "Taiwan"
## [70] "The former Yugoslav Republic of Macedonia"
## [71] "Tokelau"
## [72] "Turks and Caicos Islands"
## [73] "Tuvalu"
## [74] "United Republic of Tanzania"
## [75] "United States"
## [76] "United States Minor Outlying Islands"
## [77] "United States Virgin Islands"
## [78] "Wallis and Futuna Islands"
## [79] "Western Sahara"
## [80] "Yemen"
sort(setdiff(epi$country, world.spdf$NAME))
## [1] "Cabo Verde" "Dem. Rep. Congo"
## [3] "Eswatini" "Iran"
## [5] "Laos" "Micronesia"
## [7] "Moldova" "Myanmar"
## [9] "North Macedonia" "Republic of Congo"
## [11] "South Korea" "Tanzania"
## [13] "United States of America"
Now, collate the naming schemes between the epi data.frame and the world.spdf SpatialPolygonsDataFrame with the helper function collateCountryNames().
epi.world.spdf <- collateCountryNames(epi, world.spdf)
list2env(epi.world.spdf, .GlobalEnv)
## <environment: R_GlobalEnv>
rm(epi.world.spdf)
Lastly, control if all EPI countries are present in the world map…
sort(setdiff(epi$country, world.spdf$NAME))
## character(0)
… as the set difference is now equal to zero (i.e., all items in epi$country are present also in world.spdf$NAME).
Countries without an EPI are miniature states, islands or states currently in an armed conflict (some of the island belong to a nation; e.g., Greenland to Denkmark and Svalbard to Norway. For now, keep it as is and later, these islands could possibly be merged into their mainlands).
sort(setdiff(world.spdf$NAME, epi$country))
## [1] "Åland Islands"
## [2] "American Samoa"
## [3] "Andorra"
## [4] "Anguilla"
## [5] "Antarctica"
## [6] "Aruba"
## [7] "Bermuda"
## [8] "Bouvet Island"
## [9] "British Indian Ocean Territory"
## [10] "British Virgin Islands"
## [11] "Cayman Islands"
## [12] "Christmas Island"
## [13] "Cocos (Keeling) Islands"
## [14] "Cook Islands"
## [15] "Falkland Islands (Malvinas)"
## [16] "Faroe Islands"
## [17] "French Guiana"
## [18] "French Polynesia"
## [19] "French Southern and Antarctic Lands"
## [20] "Gibraltar"
## [21] "Greenland"
## [22] "Guadeloupe"
## [23] "Guam"
## [24] "Guernsey"
## [25] "Heard Island and McDonald Islands"
## [26] "Holy See (Vatican City)"
## [27] "Hong Kong"
## [28] "Isle of Man"
## [29] "Jersey"
## [30] "Libyan Arab Jamahiriya"
## [31] "Liechtenstein"
## [32] "Macau"
## [33] "Martinique"
## [34] "Mayotte"
## [35] "Monaco"
## [36] "Montserrat"
## [37] "Nauru"
## [38] "Netherlands Antilles"
## [39] "New Caledonia"
## [40] "Niue"
## [41] "Norfolk Island"
## [42] "North Korea"
## [43] "Northern Mariana Islands"
## [44] "Palau"
## [45] "Palestine"
## [46] "Pitcairn Islands"
## [47] "Puerto Rico"
## [48] "Reunion"
## [49] "Saint Barthelemy"
## [50] "Saint Helena"
## [51] "Saint Kitts and Nevis"
## [52] "Saint Martin"
## [53] "Saint Pierre and Miquelon"
## [54] "San Marino"
## [55] "Somalia"
## [56] "South Georgia South Sandwich Islands"
## [57] "Svalbard"
## [58] "Syrian Arab Republic"
## [59] "Taiwan"
## [60] "Tokelau"
## [61] "Turks and Caicos Islands"
## [62] "Tuvalu"
## [63] "United States Minor Outlying Islands"
## [64] "United States Virgin Islands"
## [65] "Wallis and Futuna Islands"
## [66] "Western Sahara"
## [67] "Yemen"
world.spdf <- merge(world.spdf,
epi[, (colnames(epi) %in% c("country", "EPI.new", "region",
"X2019"))],
by.x = "NAME", by.y = "country")
names(world.spdf)
## [1] "NAME" "FIPS" "ISO2" "ISO3" "UN" "AREA"
## [7] "POP2005" "REGION" "SUBREGION" "LON" "LAT" "EPI.new"
## [13] "region" "X2019"
world.spdf@data$POP2005[which(world.spdf@data$POP2005 == 0)] <- NA
world.spdf@data$POP2005 <- as.integer(world.spdf@data$POP2005) / 1000000 %>%
round(2)
The EPI color scheme:
epi.cols <- colorNumeric(palette = "viridis", domain = world.spdf@data$EPI.new,
na.color = "gray", reverse = TRUE)
The region color scheme:
region.cols <- colorFactor(palette = "plasma", domain = world.spdf@data$region,
na.color = "gray")
tooltips <- paste("Country:", world.spdf@data$NAME, "<br/>",
"Region:", world.spdf@data$region, "<br/r>",
# "Area:", world.spdf@data$AREA, "<br/>",
"Population:", round(world.spdf@data$POP2005, 2), "M",
"<br/>",
"EPI:", world.spdf@data$EPI.new) %>%
lapply(htmltools::HTML)
cp.map <- leaflet(world.spdf) %>%
addTiles() %>%
setView(lat = 20, lng = 0, zoom = 2.5) %>%
addPolygons(fillColor = ~ epi.cols(EPI.new), stroke = TRUE, group = "EPI",
fillOpacity = 0.9, color = "white", weight = 0.5,
label = tooltips,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "13px",
direction = "auto")) %>%
addPolygons(fillColor = ~ region.cols(region), stroke = TRUE,
group = "Region", fillOpacity = 0.5, color = "white",
weight = 0.5, label = tooltips,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "13px",
direction = "auto")) %>%
addLegend(pal = epi.cols, values = ~ EPI.new, opacity = 0.9,
title = "EPI", position = "topleft", group = "EPI") %>%
addLegend(pal = region.cols, values = ~ region, opacity = 0.5,
title = "", position = "topleft", group = "Region") %>%
addLayersControl(overlayGroups = c("EPI", "Region"),
options = layersControlOptions(collapsed = FALSE)) %>%
hideGroup("Region")
cp.map